Project Viva: Cord blood DNAm and metals EWAS

Required Packages and functions

# libraries
library(sas7bdat)
library(sva)
library(reshape)
library(ggplot2)
library(gridExtra)
library(knitr)
library(dplyr)
library(kableExtra)
library(tidyverse)
library(minfi)
library(stringr)
library(limma)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
library(minfi)
library(DMRcate)
library(UpSetR)
library(reshape)
library(corrplot)
library(factoextra)
library(ENmix)
library(missMethyl)

# loading annotation
anno = getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
anno = data.frame(anno)

# loading functions
setwd('/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local')
source('EWAS functions.R')

Load DNAm and pheno data

load("/Users/annebozack/Box/NIEHS-R01 ONES/Methylation Data/CordBlood_ComBat_Betas_Mvals_filteredProbes_metalAnalysis.RData")

dim(pDatcordMetal)
# 361 164

dim(ComBat.Mvalues.Metals)
# 394460    361

rownames(pDatcordMetal) = pDatcordMetal$samplename

EWAS

As

DMP_As_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'As_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  11440 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  19530 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.065045 
 
# Number of DMRs identified:   1 

# rm(DMP_As_cord); gc()

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$As_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.1888243

TestCDFs2(Z = Z, y = y, X = X)
# 0.4438389

# rm(DMP_As_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
12	9217389	9217859	0	470	0	0	9

No adjustment for fish consumption

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As/As_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As/As_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As/As_log2_manhattan_DMP_adj.png")

As, female

pDatcordMetal_F = pDatcordMetal[pDatcordMetal$female_d == 1,]
pDatcordMetal_F$race_child2 = as.factor(as.numeric(pDatcordMetal_F$race_child2))

DMP_As_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'As_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  17949 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  16668 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.9807579 

# Number of DMRs identified:   2 

# Global DNAm
ComBat.Betas.Metals_F = ComBat.Betas.Metals[,colnames(ComBat.Betas.Metals) %in% rownames(pDatcordMetal_F)]
ComBat.Betas.Metals_F = ComBat.Betas.Metals_F[,match(rownames(pDatcordMetal_F), colnames(ComBat.Betas.Metals_F))]

Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$As_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.2559703

TestCDFs2(Z = Z, y = y, X = X)
# 0.6246159

# rm(DMP_As_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
3	186490655	186490915	0	260	0	1.8e-06	5
11	34460106	34460386	0	280	0	4.7e-06	7

No adjustment for fish consumption

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_F/As_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_F/As_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_F/As_log2_manhattan_DMP_adj.png")

As, male

pDatcordMetal_M = pDatcordMetal[pDatcordMetal$female_d == 0,]
pDatcordMetal_M$race_child2 = as.factor(as.numeric(pDatcordMetal_M$race_child2))

DMP_As_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'As_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  12218 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  17743 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.011027 
 
# Number of DMRs identified:   6 

# Global DNAm
ComBat.Betas.Metals_M = ComBat.Betas.Metals[,colnames(ComBat.Betas.Metals) %in% rownames(pDatcordMetal_M)]
ComBat.Betas.Metals_M = ComBat.Betas.Metals_M[,match(rownames(pDatcordMetal_M), colnames(ComBat.Betas.Metals_M))]

Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$As_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.5073422

TestCDFs2(Z = Z, y = y, X = X)
# 0.4807815

# rm(DMP_As_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
12	9217389	9217859	0.0000000	470	0.0000000	0.0000000	9
20	3051953	3052345	0.0000000	392	0.0000000	0.0000001	9
11	18433499	18433683	0.0000000	184	0.0000001	0.0000770	4
11	2890628	2890670	0.0000178	42	0.0000214	0.1539887	6
11	1036676	1036766	0.0000179	90	0.0000214	0.0753032	3
6	32063990	32064032	0.0060759	42	0.0060759	1.0000000	2

No adjustment for fish consumption

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_M/As_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_M/As_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/As_M/As_log2_manhattan_DMP_adj.png")

Ba

DMP_Ba_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Ba_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  21408 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  18967 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.006289 
 
# Number of identified DMR:  0

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Ba_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.9650461

# rm(DMP_Ba_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
NA	NA	NA	NA	NA	NA	NA	NA

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba/Ba_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba/Ba_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba/Ba_log2_manhattan_DMP_adj.png")

Ba, female

DMP_Ba_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Ba_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  16994 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  13428 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.7990795 

# Number of DMRs identified:   2 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Ba_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.864574

TestCDFs2(Z = Z, y = y, X = X)
# 0.9088558

# rm(DMP_Ba_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
10	135051255	135051581	0	326	0	2e-06	7
8	144635259	144635610	0	351	0	5e-06	9

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_F/Ba_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_F/Ba_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_F/Ba_log2_manhattan_DMP_adj.png")

Ba, male

DMP_Ba_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Ba_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  16173 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  14282 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.865515 
 
# Number of DMRs identified:   2 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Ba_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.6927273

TestCDFs2(Z = Z, y = y, X = X)
# 0.8277251

# rm(DMP_Ba_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
4	169239618	169240010	0	392	0	0e+00	7
10	135278716	135279147	0	431	0	2e-07	5

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_M/Ba_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_M/Ba_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Ba_M/Ba_log2_manhattan_DMP_adj.png")

Cd

DMP_Cd_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Cd_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  13284 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  16155 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.9008795 

# Number of DMRs identified:   5 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Cd_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.9348243

# rm(DMP_Cd_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
3	156323951	156324118	0.0000000	167	0.0000000	0.0000030	3
13	110319561	110319607	0.0000000	46	0.0000000	0.0000233	3
17	46685291	46685448	0.0000000	157	0.0000000	0.0000207	5
2	200468625	200468832	0.0000001	207	0.0000001	0.0001106	3
6	30039141	30039175	0.0095070	34	0.0095070	1.0000000	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd/Cd_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd/Cd_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd/Cd_log2_manhattan_DMP_adj.png")

Cd, female

DMP_Cd_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Cd_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  19177 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  18135 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.9732893 
 
# Number of DMRs identified:   2 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Cd_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.6678424

TestCDFs2(Z = Z, y = y, X = X)
# 0.7811745

# rm(DMP_Cd_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
7	27225810	27225897	2.0e-07	87	4.0e-07	0.0009656	4
2	200468625	200468728	4.5e-06	103	4.5e-06	0.0171663	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_F/Cd_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_F/Cd_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_F/Cd_log2_manhattan_DMP_adj.png")

Cd, male

DMP_Cd_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Cd_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  8842 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  15768 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.9683207 
 
# Number of DMRs identified:   1

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Cd_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.4173395
 
# rm(DMP_Cd_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
19	37742738	37742956	0	218	0	1e-07	4

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_M/Cd_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_M/Cd_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cd_M/Cd_log2_manhattan_DMP_adj.png")

Cr

DMP_Cr_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Cr_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  19469 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  19356 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.9710716 

# Number of DMRs identified:   5 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Cr_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.8160506

TestCDFs2(Z = Z, y = y, X = X)
# 0.8524831

# rm(DMP_Cr_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
8	13373032	13373141	0.0e+00	109	0.0e+00	0.0000063	3
6	33245700	33246008	0.0e+00	308	0.0e+00	0.0000023	13
10	114713023	114713187	0.0e+00	164	0.0e+00	0.0000259	3
17	46681110	46681401	0.0e+00	291	1.0e-07	0.0000623	6
17	46641862	46642011	8.8e-05	149	8.8e-05	0.2079115	2

No adjustment for fish consumption

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr/Cr_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr/Cr_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr/Cr_log2_manhattan_DMP_adj.png")

Cr, female

DMP_Cr_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Cr_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  15600 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  13645 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.8363976 
 
# Number of DMRs identified:   2 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Cr_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.9279794

# rm(DMP_Cr_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
8	13373032	13373141	0e+00	109	0e+00	0.0000145	3
6	31148369	31148666	1e-07	297	1e-07	0.0001264	12

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_F/Cr_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_F/Cr_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_F/Cr_log2_manhattan_DMP_adj.png")

Cr, male

DMP_Cr_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Cr_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  13268 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  15910 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.909677 

# Number of identified DMR:  0

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Cr_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.5576744

TestCDFs2(Z = Z, y = y, X = X)
# 0.7528378

rm(DMP_Cr_M_cord); gc()

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_M/Cr_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_M/Cr_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cr_M/Cr_log2_manhattan_DMP_adj.png")

Cs

DMP_Cs_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Cs_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  28868 
# Unadjusted, FDR<0.05:  5 
# Unadjusted, pBonf<0.05:  3 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  19239 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.018459 

# Number of DMRs identified:   18 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Cs_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.5068162

TestCDFs2(Z = Z, y = y, X = X)
# 0.5042182

# rm(DMP_Cs_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
4	174421376	174422626	0.0000000	1250	0.0000000	0.0000000	7
6	33280051	33280478	0.0000000	427	0.0000000	0.0000000	14
6	28601268	28601519	0.0000000	251	0.0000000	0.0000001	11
6	32847440	32847845	0.0000000	405	0.0000000	0.0000003	17
6	30881315	30881842	0.0000000	527	0.0000000	0.0000002	21
17	17603530	17603837	0.0000000	307	0.0000000	0.0000008	4
20	57427411	57427762	0.0000000	351	0.0000000	0.0000016	15
1	75590911	75591029	0.0000000	118	0.0000000	0.0000167	3
12	44152508	44152940	0.0000000	432	0.0000000	0.0000091	11
12	1025528	1025755	0.0000001	227	0.0000001	0.0000898	3
12	54763210	54763433	0.0000001	223	0.0000001	0.0000934	3
19	10736005	10736117	0.0000001	112	0.0000001	0.0002107	5
3	141087186	141087363	0.0000001	177	0.0000001	0.0001661	5
8	43131259	43131656	0.0000001	397	0.0000001	0.0000851	5
15	69325270	69325560	0.0000001	290	0.0000001	0.0001177	5
5	78985424	78985592	0.0000002	168	0.0000002	0.0003553	9
17	46388389	46388465	0.0000002	76	0.0000003	0.0012942	3
6	31651019	31651029	0.0061778	10	0.0061778	1.0000000	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs/Cs_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs/Cs_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs/Cs_log2_manhattan_DMP_adj.png")

Cs, female

DMP_Cs_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Cs_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  14228 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  11077 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.7878935 
 
# Number of DMRs identified:   4 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Cs_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.683642

TestCDFs2(Z = Z, y = y, X = X)
# 0.9583828

# rm(DMP_Cs_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
12	44152508	44152940	0.00e+00	432	0.00e+00	0.0000000	11
6	117923794	117924070	0.00e+00	276	0.00e+00	0.0000220	8
1	26233331	26233623	0.00e+00	292	1.00e-07	0.0000511	10
6	32847761	32847830	1.61e-05	69	1.61e-05	0.0879654	5

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_F/Cs_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_F/Cs_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_F/Cs_log2_manhattan_DMP_adj.png")

Cs, male

DMP_Cs_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Cs_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  23725 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  21075 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.092447 

# Number of DMRs identified:   14 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Cs_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.9380089

TestCDFs2(Z = Z, y = y, X = X)
# 0.5680556

# rm(DMP_Cs_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
6	31650734	31651362	0e+00	628	0e+00	0.0000000	21
8	43131259	43132507	0e+00	1248	0e+00	0.0000000	8
15	69325270	69325560	0e+00	290	0e+00	0.0000000	5
17	17603530	17604146	0e+00	616	0e+00	0.0000000	5
12	108634146	108634275	0e+00	129	0e+00	0.0000017	3
6	28601268	28601519	0e+00	251	0e+00	0.0000024	11
3	141087186	141087363	0e+00	177	0e+00	0.0000098	5
10	134150488	134150760	0e+00	272	1e-07	0.0000454	7
6	30881463	30881766	1e-07	303	1e-07	0.0001063	15
19	51774376	51774666	1e-07	290	1e-07	0.0001355	4
14	21148813	21148957	1e-07	144	1e-07	0.0003044	2
6	28446839	28447115	1e-07	276	2e-07	0.0002057	4
4	187422004	187422119	2e-07	115	2e-07	0.0006142	5
13	23309929	23310203	4e-07	274	4e-07	0.0005223	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_M/Cs_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_M/Cs_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cs_M/Cs_log2_manhattan_DMP_adj.png")

Cu

DMP_Cu_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Cu_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  31830 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  57348 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.845378 

# Number of DMRs identified:   4 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Cu_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.3984576

TestCDFs2(Z = Z, y = y, X = X)
# 0.05580438

# rm(DMP_Cu_cord); gc()

FDR-significant DMPs

x
NA

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
4	174421376	174422626	0.0000000	1250	0.0000000	0.0000000	7
6	33280051	33280478	0.0000000	427	0.0000000	0.0000000	14
6	28601268	28601519	0.0000000	251	0.0000000	0.0000001	11
6	32847440	32847845	0.0000000	405	0.0000000	0.0000003	17
6	30881315	30881842	0.0000000	527	0.0000000	0.0000002	21
17	17603530	17603837	0.0000000	307	0.0000000	0.0000008	4
20	57427411	57427762	0.0000000	351	0.0000000	0.0000016	15
1	75590911	75591029	0.0000000	118	0.0000000	0.0000167	3
12	44152508	44152940	0.0000000	432	0.0000000	0.0000091	11
12	1025528	1025755	0.0000001	227	0.0000001	0.0000898	3
12	54763210	54763433	0.0000001	223	0.0000001	0.0000934	3
19	10736005	10736117	0.0000001	112	0.0000001	0.0002107	5
3	141087186	141087363	0.0000001	177	0.0000001	0.0001661	5
8	43131259	43131656	0.0000001	397	0.0000001	0.0000851	5
15	69325270	69325560	0.0000001	290	0.0000001	0.0001177	5
5	78985424	78985592	0.0000002	168	0.0000002	0.0003553	9
17	46388389	46388465	0.0000002	76	0.0000003	0.0012942	3
6	31651019	31651029	0.0061778	10	0.0061778	1.0000000	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu/Cu_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu/Cu_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu/Cu_log2_manhattan_DMP_adj.png")

Cu, female

DMP_Cu_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Cu_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  34435 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  17137 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.061581 
 
# Number of DMRs identified:   3 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Cu_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.2253006

TestCDFs2(Z = Z, y = y, X = X)
# 0.2135721

# rm(DMP_Cu_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
20	57427442	57427830	0.0000000	388	0.0000000	0.0000017	15
7	1250087	1250182	0.0000576	95	0.0000864	0.2126146	3
11	368637	368712	0.0067564	75	0.0067564	1.0000000	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_F/Cu_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_F/Cu_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_F/Cu_log2_manhattan_DMP_adj.png")

Cu, male

DMP_Cu_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Cu_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  10164 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  57164 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.911863 
 
# Number of DMRs identified:   6 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Cu_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.2307306

# rm(DMP_Cu_M); gc()

FDR-significant DMPs

x
NA

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
5	83016999	83017184	0.0000000	185	0.0000000	0.0000001	4
3	46759437	46759698	0.0000000	261	0.0000000	0.0000021	7
1	108023248	108023482	0.0000001	234	0.0000001	0.0001015	5
12	54673866	54674009	0.0000036	143	0.0000053	0.0097873	4
6	32847761	32847830	0.0000643	69	0.0000772	0.3076772	5
7	94953876	94953956	0.0049862	80	0.0049862	1.0000000	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_M/Cu_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_M/Cu_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Cu_M/Cu_log2_manhattan_DMP_adj.png")

Hg

DMP_Hg_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Hg_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg')

# Unadjusted, N =  358 
# Unadjusted, p<0.05:  38119 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  358 
# Adjusted, p<0.05:  13462 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.7929906 
 
# Number of DMRs identified:   2 

# Global DNAm
pDatcordMetal_Hg = pDatcordMetal[!is.na(pDatcordMetal$Hg_log2),]
dim(pDatcordMetal_Hg)
# 358 164

ComBat.Betas.Metals_Hg = ComBat.Betas.Metals[,colnames(ComBat.Betas.Metals) %in% rownames(pDatcordMetal_Hg)]
ComBat.Betas.Metals_Hg = ComBat.Betas.Metals_Hg[,match(rownames(pDatcordMetal_Hg), colnames(ComBat.Betas.Metals_Hg))]

Z = as.matrix(ComBat.Betas.Metals_Hg)
y = pDatcordMetal_Hg$Hg_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_Hg)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.845803

# rm(DMP_Hg_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
8	143859668	143859990	0	322	0	0.00e+00	7
6	30039373	30039548	0	175	0	1.99e-05	12

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg/Hg_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg/Hg_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg/Hg_log2_manhattan_DMP_adj.png")

Hg, female

DMP_Hg_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Hg_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_F')

# Unadjusted, N =  167 
# Unadjusted, p<0.05:  11528 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  167 
# Adjusted, p<0.05:  10298 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.7271097 
 
# Number of DMRs identified:   1 

# Global DNAm
pDatcordMetal_HgF = pDatcordMetal_Hg[pDatcordMetal_Hg$female_d == 1,]
dim(pDatcordMetal_HgF)
# 167 164

ComBat.Betas.Metals_HgF = ComBat.Betas.Metals[,colnames(ComBat.Betas.Metals) %in% rownames(pDatcordMetal_HgF)]
ComBat.Betas.Metals_HgF = ComBat.Betas.Metals_HgF[,match(rownames(pDatcordMetal_HgF), colnames(ComBat.Betas.Metals_HgF))]

Z = as.matrix(ComBat.Betas.Metals_HgF)
y = pDatcordMetal_HgF$Hg_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_HgF)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.9751191

# rm(DMP_Hg_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
6	31508105	31508137	3.95e-05	32	3.95e-05	0.3855435	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_F/Hg_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_F/Hg_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_F/Hg_log2_manhattan_DMP_adj.png")

Hg, male

DMP_Hg_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Hg_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_M')

# Unadjusted, N =  191 
# Unadjusted, p<0.05:  32177 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  191 
# Adjusted, p<0.05:  12873 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.8093155 
 
# Number of DMRs identified:   3 

# Global DNAm
pDatcordMetal_HgM = pDatcordMetal_Hg[pDatcordMetal_Hg$female_d == 0,]
dim(pDatcordMetal_HgM)
# 191 164

ComBat.Betas.Metals_HgM = ComBat.Betas.Metals[,colnames(ComBat.Betas.Metals) %in% rownames(pDatcordMetal_HgM)]
ComBat.Betas.Metals_HgM = ComBat.Betas.Metals_HgM[,match(rownames(pDatcordMetal_HgM), colnames(ComBat.Betas.Metals_HgM))]

Z = as.matrix(ComBat.Betas.Metals_HgM)
y = pDatcordMetal_HgM$Hg_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_HgM)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.8570357

# rm(DMP_Hg_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
8	143859668	143859990	0	322	0	0e+00	7
6	31650785	31651362	0	577	0	2e-07	19
6	30039141	30039548	0	407	0	4e-07	15

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_M/Hg_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_M/Hg_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Hg_M/Hg_log2_manhattan_DMP_adj.png")

Mg

DMP_Mg_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Mg_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  20828 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  19876 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.067098 

# Number of DMRs identified:   5 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Mg_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.469153

TestCDFs2(Z = Z, y = y, X = X)
# 0.3509856

# rm(DMP_Mg_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
20	57427169	57427973	0.0e+00	804	0.0e+00	0.0000000	24
6	30039174	30039548	0.0e+00	374	0.0e+00	0.0000000	13
6	32063990	32064258	0.0e+00	268	1.0e-07	0.0000461	12
1	2980037	2980163	2.5e-06	126	3.1e-06	0.0077603	2
11	368564	368712	6.6e-06	148	6.6e-06	0.0175234	7

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg/Mg_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg/Mg_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg/Mg_log2_manhattan_DMP_adj.png")

Mg, female

DMP_Mg_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Mg_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  21465 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  15005 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.860424 
 
# Number of DMRs identified:   5

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Mg_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.4382674

TestCDFs2(Z = Z, y = y, X = X)
# 0.5358737

# rm(DMP_Mg_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
17	48912264	48912621	0.00e+00	357	0.00e+00	0.0000000	9
20	57427442	57427762	0.00e+00	320	1.00e-07	0.0000298	13
2	121223739	121224009	1.00e-07	270	1.00e-07	0.0001158	6
14	105944941	105945022	1.60e-06	81	2.00e-06	0.0076408	2
19	45976119	45976195	4.01e-05	76	4.01e-05	0.1879917	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_F/Mg_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_F/Mg_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_F/Mg_log2_manhattan_DMP_adj.png")

Mg, male

DMP_Mg_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Mg_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  38543 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  40446 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.605054 
 
# Number of DMRs identified:   2 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Mg_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.6190037

TestCDFs2(Z = Z, y = y, X = X)
# 0.3666586

# rm(DMP_Mg_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
19	57742111	57742423	0.0000000	312	0.0000000	0.0000000	7
22	22901829	22901880	0.0008779	51	0.0008779	0.9988789	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_M/Mg_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_M/Mg_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mg_M/Mg_log2_manhattan_DMP_adj.png")

Mn

DMP_Mn_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Mn_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  9732 
# Unadjusted, FDR<0.05:  1 
# Unadjusted, pBonf<0.05:  1 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  10167 
# Adjusted, FDR<0.05:  1 
# Adjusted, pBonf<0.05:  1 
# Adjusted, lambda:  0.7325852 
 
# Number of DMRs identified:   2 

# Pathway analysis

GO_Mn = gometh(sig.cpg = DMP_Mn_cord[[3]]$cpg[DMP_Mn_cord[[3]]$adj.P.Val < 0.05], all.cpg = DMP_Mn_cord[[3]]$cpg, collection = "GO")
table(GO_Mn$P.DE < 0.05)
# FALSE  TRUE 
# 22583    27 
table(GO_Mn$FDR < 0.05)
# FALSE 
# 22610 

KEGG_Mn = gometh(sig.cpg = DMP_Mn_cord[[3]]$cpg[DMP_Mn_cord[[3]]$adj.P.Val < 0.05], all.cpg = DMP_Mn_cord[[3]]$cpg, collection = "KEGG")
table(KEGG_Mn$P.DE < 0.05)
# FALSE 
  # 342 
  
# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Mn_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.9723614

# rm(DMP_Mn_cord); gc()

FDR-significant DMPs

	logFC_CI	AveExpr	P.Value	adj.P.Val	adj.P.Val.bonf	chr	pos	gene
cg02042823	0.43 (0.31 ,0.55)	5.75	0	9.3e-06	0	chr16	6714429	A2BP1;A2BP1

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
6	32063725	32064161	0.0000000	436	0.0000000	0.0000038	13
7	1250125	1250182	0.0017739	57	0.0017739	0.9999954	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn/Mn_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn/Mn_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn/Mn_log2_manhattan_DMP_adj.png")

Mn, female

DMP_Mn_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Mn_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  26444 
# Unadjusted, FDR<0.05:  3 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  13233 
# Adjusted, FDR<0.05:  9 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.8939725 
 
# Number of DMRs identified:   7 

GO_Mn_F = gometh(sig.cpg = DMP_Mn_F_cord[[3]]$cpg[DMP_Mn_F_cord[[3]]$adj.P.Val < 0.05], all.cpg = DMP_Mn_F_cord[[3]]$cpg, collection = "GO")
table(GO_Mn_F$P.DE < 0.05)
# FALSE  TRUE 
# 22481   129 
table(GO_Mn_F$FDR < 0.05)
# FALSE 
# 22610 

KEGG_Mn_F = gometh(sig.cpg = DMP_Mn_F_cord[[3]]$cpg[DMP_Mn_F_cord[[3]]$adj.P.Val < 0.05], all.cpg = DMP_Mn_cord_F[[3]]$cpg, collection = "KEGG")
table(KEGG_Mn_F$P.DE < 0.05)
# FALSE  TRUE 
  # 334     8 
table(KEGG_Mn_F$FDR < 0.05)
# FALSE 
  # 342 
  
write.csv(KEGG_Mn_F, '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_F/KEGG.csv')

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Mn_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.3499748

TestCDFs2(Z = Z, y = y, X = X)
# 0.7004179

# rm(DMP_Mn_F_cord); gc()

FDR-significant DMPs

	logFC_CI	AveExpr	P.Value	adj.P.Val	adj.P.Val.bonf	chr	pos	gene
cg00954161	0.42 (0.27 ,0.57)	5.99	2e-07	0.0341269	0.0801	chr1	3696925	LRRC47
cg01744822	-0.22 (-0.31 ,-0.14)	-4.36	9e-07	0.0442988	0.3421	chr16	73100510
cg08904630	0.21 (0.13 ,0.29)	5.33	9e-07	0.0442988	0.3514	chr10	71490427
cg11161853	-0.14 (-0.19 ,-0.08)	-5.29	9e-07	0.0442988	0.3738	chr3	67705044	SUCLG2
cg15712310	-0.19 (-0.26 ,-0.12)	-1.62	2e-07	0.0341269	0.0717	chr16	73100790
cg19908812	-0.17 (-0.24 ,-0.1)	-6.02	9e-07	0.0442988	0.3659	chr4	164253006	NPY1R
cg22799518	0.52 (0.33 ,0.71)	6.34	3e-07	0.0341269	0.1365	chr12	56988862	RBMS2
cg23903787	0.34 (0.22 ,0.47)	3.02	3e-07	0.0341269	0.1141	chr4	3527371	LRPAP1
cg26462130	0.39 (0.24 ,0.53)	6.11	1e-06	0.0442988	0.3987	chr7	2052961	MAD1L1;MAD1L1;MAD1L1

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
16	73100425	73100790	0e+00	365	0e+00	0.0000000	3
3	148804271	148804556	0e+00	285	0e+00	0.0000004	8
3	67704889	67705285	0e+00	396	0e+00	0.0000005	6
6	161561030	161561121	0e+00	91	0e+00	0.0000071	3
7	3019158	3019382	0e+00	224	0e+00	0.0000054	4
7	27170716	27171051	0e+00	335	0e+00	0.0000047	9
13	97646638	97646765	5e-07	127	5e-07	0.0016533	4

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_F/Mn_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_F/Mn_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_F/Mn_log2_manhattan_DMP_adj.png")

Mn, male

DMP_Mn_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Mn_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  9543 
# Unadjusted, FDR<0.05:  2 
# Unadjusted, pBonf<0.05:  2 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  12799 
# Adjusted, FDR<0.05:  2 
# Adjusted, pBonf<0.05:  2 
# Adjusted, lambda:  0.8388281 
 
# Number of DMRs identified:   4

GO_Mn_M = gometh(sig.cpg = DMP_Mn_M_cord[[3]]$cpg[DMP_Mn_M_cord[[3]]$adj.P.Val < 0.05], all.cpg = DMP_Mn_M_cord[[3]]$cpg, collection = "GO")
table(GO_Mn_M$P.DE < 0.05)
# FALSE  TRUE 
# 22586    24 
table(GO_Mn_M$FDR < 0.05)
# FALSE 
# 22610 

KEGG_Mn_M = gometh(sig.cpg = DMP_Mn_M_cord[[3]]$cpg[DMP_Mn_M_cord[[3]]$adj.P.Val < 0.05], all.cpg = DMP_Mn_M_cord[[3]]$cpg, collection = "KEGG")
table(KEGG_Mn_M$P.DE < 0.05)
# FALSE 
  # 342 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Mn_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.3561523

TestCDFs2(Z = Z, y = y, X = X)
# 0.5219932

# rm(DMP_Mn_M_cord); gc()

FDR-significant DMPs

	logFC_CI	AveExpr	P.Value	adj.P.Val	adj.P.Val.bonf	chr	pos	gene
cg02042823	0.51 (0.36 ,0.66)	5.71	0	0.0001241	0.0001	chr16	6714429	A2BP1;A2BP1
cg03763518	-0.29 (-0.39 ,-0.19)	-3.51	0	0.0091811	0.0184	chr1	150245044	C1orf54

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
20	36148603	36149271	0e+00	668	0e+00	0.0000000	31
7	1250087	1250273	0e+00	186	0e+00	0.0000000	7
15	99789621	99789855	0e+00	234	0e+00	0.0000030	5
1	75198817	75199117	1e-07	300	1e-07	0.0001093	6

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_M/Mn_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_M/Mn_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Mn_M/Mn_log2_manhattan_DMP_adj.png")

Pb

DMP_Pb_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Pb_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  25856 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  14691 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.8763417 
 
# Number of identified DMR:  0

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Pb_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.2242514

TestCDFs2(Z = Z, y = y, X = X)
# 0.9723614

# rm(DMP_Pb_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
NA	NA	NA	NA	NA	NA	NA	NA

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb/Pb_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb/Pb_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb/Pb_log2_manhattan_DMP_adj.png")

Pb, female

DMP_Pb_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Pb_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  16597 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  20414 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.132739 
 
# Number of identified DMR:  0

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Pb_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.07160096

TestCDFs2(Z = Z, y = y, X = X)
# 0.2744584

# rm(DMP_Pb_F_cord); gc()

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_F/Pb_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_F/Pb_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_F/Pb_log2_manhattan_DMP_adj.png")

Pb, male

DMP_Pb_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Pb_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  19998 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  13522 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.8591965 
 
# Number of DMRs identified:   2

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Pb_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 1

TestCDFs2(Z = Z, y = y, X = X)
# 0.9140094 

# rm(DMP_Pb_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
14	100203941	100204258	0	317	0	1.9e-06	6
3	122640777	122641144	0	367	0	6.5e-06	4

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_M/Pb_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_M/Pb_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Pb_M/Pb_log2_manhattan_DMP_adj.png")

Se

DMP_Se_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Se_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  23447 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  15897 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.8643957 

# Number of DMRs identified:   2 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Se_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.1214307

TestCDFs2(Z = Z, y = y, X = X)
# 0.2808225

# rm(DMP_Se_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
8	144635259	144635610	0	351	0	7.90e-06	9
4	74847645	74847829	0	184	0	4.18e-05	7

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se/Se_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se/Se_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se/Se_log2_manhattan_DMP_adj.png")

Se, female

DMP_Se_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Se_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_F')

# # Unadjusted, N =  169 
# Unadjusted, p<0.05:  11702 
# Unadjusted, FDR<0.05:  1 
# Unadjusted, pBonf<0.05:  1 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  10751 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  0.717438 
 
# Number of DMRs identified:   4 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Se_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.6522482

TestCDFs2(Z = Z, y = y, X = X)
# 0.876843

# rm(DMP_Se_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
8	144635259	144636113	0.0000000	854	0.0000000	0.0000000	12
15	75019069	75019376	0.0000000	307	0.0000000	0.0000318	10
19	1510493	1510692	0.0000000	199	0.0000000	0.0000645	4
18	23713837	23714009	0.0002329	172	0.0002329	0.4138549	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_F/Se_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_F/Se_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_F/Se_log2_manhattan_DMP_adj.png")

Se, male

DMP_Se_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Se_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  28127 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  23097 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.035939 
 
# Number of DMRs identified:   1 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Se_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.04829198

TestCDFs2(Z = Z, y = y, X = X)
# 0.1484135 

# rm(DMP_Se_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
19	57742254	57742423	0	169	0	2.7e-06	6

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_M/Se_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_M/Se_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Se_M/Se_log2_manhattan_DMP_adj.png")

Zn

DMP_Zn_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal, var = 'Zn_log2', covar = c('female_d', 'race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn')

# Unadjusted, N =  361 
# Unadjusted, p<0.05:  25089 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  361 
# Adjusted, p<0.05:  37705 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.490786 
 
# Number of DMRs identified:   6 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals)
y = pDatcordMetal$Zn_log2
X = model.matrix(~ female_d + race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.6606383

TestCDFs2(Z = Z, y = y, X = X)
# 0.2234394

# rm(DMP_Zn); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
3	48632483	48632783	0.00e+00	300	0.00e+00	0.0000004	8
17	79503641	79503877	0.00e+00	236	0.00e+00	0.0000010	4
12	1973871	1974168	0.00e+00	297	0.00e+00	0.0000053	3
11	92702372	92702912	0.00e+00	540	0.00e+00	0.0000053	8
20	25129506	25129562	5.00e-06	56	6.00e-06	0.0348380	5
6	32165175	32165200	7.09e-05	25	7.09e-05	0.6735463	3

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn/Zn_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn/Zn_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn/Zn_log2_manhattan_DMP_adj.png")

Zn, female

DMP_Zn_F_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_F, var = 'Zn_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_F')

# Unadjusted, N =  169 
# Unadjusted, p<0.05:  35062 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  169 
# Adjusted, p<0.05:  17286 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.062488 
 
# Number of DMRs identified:   1 

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_F)
y = pDatcordMetal_F$Zn_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_F)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.5894776

TestCDFs2(Z = Z, y = y, X = X)
# 0.2311608

# rm(DMP_Zn_F_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
3	138763669	138763910	0	241	0	0	5

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_F/Zn_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_F/Zn_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_F/Zn_log2_manhattan_DMP_adj.png")

Zn, male

DMP_Zn_M_cord = run_EWAS(DNAm = ComBat.Mvalues.Metals, pheno = pDatcordMetal_M, var = 'Zn_log2', covar = c('race_child2', 'gestage_wks_deliv_d', 'age_mom_enroll_d', 'bmi_mom_prepreg_d', 'coll_grad', 'nullip', 'gt70k', 'smk_preg', 'Bcell_GS_cb', 'CD4T_GS_cb', 'CD8T_GS_cb', 'Gran_GS_cb', 'Mono_GS_cb', 'NK_GS_cb', 'nRBC_GS_cb'), anno = anno, path = '/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_M')

# Unadjusted, N =  192 
# Unadjusted, p<0.05:  10405 
# Unadjusted, FDR<0.05:  0 
# Unadjusted, pBonf<0.05:  0 
 
# Adjusted, N =  192 
# Adjusted, p<0.05:  22001 
# Adjusted, FDR<0.05:  0 
# Adjusted, pBonf<0.05:  0 
# Adjusted, lambda:  1.099038 
 
# Number of DMRs identified:   4

# Global DNAm
Z = as.matrix(ComBat.Betas.Metals_M)
y = pDatcordMetal_M$Zn_log2
X = model.matrix(~ race_child2 + gestage_wks_deliv_d + age_mom_enroll_d + bmi_mom_prepreg_d + coll_grad + nullip + gt70k + smk_preg + Bcell_GS_cb + CD4T_GS_cb + CD8T_GS_cb + Gran_GS_cb + Mono_GS_cb + NK_GS_cb + nRBC_GS_cb, data = pDatcordMetal_M)[,-1]

TestDensities2(Z = Z, y = y, X = X)
# 0.8342179

TestCDFs2(Z = Z, y = y, X = X)
# 0.7287762 

# rm(DMP_Zn_M_cord); gc()

DMRs

chr	start	end	p	length	fdr	sidak	nprobe
13	110319561	110319607	0.0000000	46	0.0000000	4.29e-05	3
5	112824496	112824765	0.0000000	269	0.0000000	9.90e-06	6
1	1361654	1361729	0.0000000	75	0.0000000	4.20e-05	3
6	31651019	31651029	0.0077377	10	0.0077377	1.00e+00	2

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_M/Zn_log2_QQ_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_M/Zn_log2_volcano_DMP_adj.png")

knitr::include_graphics("/Users/annebozack/Documents/Cardenas/viva_DNAm_metals_local/CordBlood/Zn_M/Zn_log2_manhattan_DMP_adj.png")